/*
 * Functions.h
 *
 *  Created on: 31/08/2012
 *      Author: Emeson J. S. Pereira
 */

#ifndef FUNCTIONS_H_
#define FUNCTIONS_H_



#endif /* FUNCTIONS_H_ */

#include <stdio.h>

//==============Functions================================
double max(register double a, register double b);
double min(register double a, register double b);
double double_sum(register double sum, register double number);
long factorial(register long n);
double distance(register double A1,register double B1,register double A2,register double B2);
double squaredroot(const double fg);
double mod(register double x, register double y);

//=======================================================

double // Max between a and b
max(register double a, register double b)
{
	if(a > b) return a;
	else return b;
}

double // Min between a and b
min(register double a, register double b)
{
	if(a < b) return a;
	else return b;
}

double //Kahan Summation Algorithm
double_sum(register double sum, register double number)
{
	register double c = 0.0;
	register double y;
	register double t;

	y = number - c;
	t = sum + y;
	c = (t - sum) - y;
	sum = t;

	return sum;
}

long //factorial of N
factorial(register long n)
{
	if(n==1) return 1;
	else return n*factorial(n-1);
}

double //distance Euclidean between two points
distance(register double A1,register double B1,register double A2,register double B2)
{
	return squaredroot((A1-A2)*(A1-A2) + (B1-B2)*(B1-B2));
}

double //Babylonian Method
squaredroot( const double fg)
{
 double n = fg / 2.0;
 double lstX = 0.0;
 while(n != lstX)
 {
 lstX = n;
 n = (n + fg/n) / 2.0;
 }
 return n;
 }
/*{//Fast Square Root Method(Babylonian Method + some manipulations on IEEE 32 bit doubleing point representation)
  union
  {
    int i;
    float x;
  } u;
  u.x = x;
  u.i = (1<<29) + (u.i >> 1) - (1<<22);

  // Two Babylonian Steps (simplified from:)
  // u.x = 0.5f * (u.x + x/u.x);
  // u.x = 0.5f * (u.x + x/u.x);
  u.x =       u.x + x/u.x;
  u.x = 0.25f*u.x + x/u.x;

  return u.x;
}*/

double
mod(register double x, register double y)
{
	register double result;
	if((x > 0 && y > 0) || (x < 0 && y < 0)) result =  x - (y*(int)(x/y));
	else if(x > 0 && y < 0) result = x - ((-y)*(1 + (-(int)(x/y))));
	else if(x < 0 && y > 0) result = x + (y*(1 + (-(int)(x/y))));

	return result;
}
